

global ROOT ""
global data "$ROOT/Data"
global dofiles "$ROOT/DoFiles"
global tables "$ROOT/Tables"
global figures "$ROOT/Figures"

use "$data\analysisfinal_adulthealth.dta", clear


tab treatment_o
*Specify control group as the base level 
fvset base 4 treatment_o



***************************************************************************
*Descriptives
***************************************************************************
table wave, c(mean diarrhoea14 mean fever14 mean cough14  mean chills14 mean vomiting14) by(male)
table wave, c(mean carry_easy mean walk_easy mean adl14_a mean adl14_n mean symptom_any ) by(male)


*******************************************************************
* Writing a program to save results from summarize command in e()
*******************************************************************
cap program drop banyy
program define banyy, eclass
args v1
sum `v1' if wg_o==0 
local mean = r(mean)
ereturn scalar mean = `mean'
end

************************************
* Anderson index of health outcomes
************************************
******************************************************
* Redefining outcomes for higher value to mean better
******************************************************
foreach y in diarrhoea14 fever14 cough14 chills14 vomiting14 adl14_a {
gen s`y' = `y'==0
replace s`y' = . if `y'==.
}

global health "sadl14_a carry_easy walk_easy sdiarrhoea14 sfever14 scough14 schills14 svomiting14"

* Doing the standardization, etc by gender

foreach y in $health {
gen std_`y' = .
qui sum `y' if male==1
local x = r(mean)
qui sum `y' if ifeed_o==0 & male==1
local s = r(sd)
replace std_`y' = (`y'-`x')/`s' if male==1

//Females now
qui sum `y' if male==0
local x = r(mean)
qui sum `y' if ifeed_o==0 & male==0
local s = r(sd)
replace std_`y' = (`y'-`x')/`s' if male==0

}

**Create index by gender
preserve
keep if male==1
* Then constructing the weighting matrix
* The following 2 lines drop missing observations
qui regress $health
keep if e(sample) 

* Next, calling mata so that we can compute the covariance matrix (which can incidentally be computed by 
* combining the standardised values above. See the formula in Appendix A of Anderson(2007).
cap mata clear
mata
st_view(X=., ., ("std_sadl14_a", "std_carry_easy", "std_walk_easy", "std_sdiarrhoea14", "std_sfever14", "std_scough14", "std_schills14", "std_svomiting14"))

Y = invsym(X'X)			/* Computing the inverse of the covariance matrix */
b = rowsum(Y)			/* Getting the sum across rows */
c= colsum(b)
/* Weights matrix is now */
weights = b/c
st_matrix("weights2", weights)
end

* Most weight is given to the diarrhoea variable, followed by vomiting and then chills
foreach y of numlist 1/8 {
local we`y' = weights2[`y',1]
}

restore

gen morb_index = `we1'*std_sadl14_a + `we2'*std_carry_easy + `we3'*std_walk_easy + `we4'*std_sdiarrhoea14 + `we5'*std_sfever14 +`we6'*std_scough14 +`we7'*std_schills14 + `we8'*std_svomiting14 if male==1

preserve
keep if male==0
* Then onstructing the weighting matrix
* The following 2 lines drop missing observations
qui regress $health
keep if e(sample) 

* Next, calling mata so that we can compute the covariance matrix (which can incidentally be computed by 
* combining the standardised values above. See the formula in Appendix A of Anderson(2007).
cap mata clear
mata
st_view(X=., ., ("std_sadl14_a", "std_carry_easy", "std_walk_easy", "std_sdiarrhoea14", "std_sfever14", "std_scough14", "std_schills14", "std_svomiting14"))

Y = invsym(X'X)			/* Computing the inverse of the covariance matrix */
b = rowsum(Y)			/* Getting the sum across rows */
c= colsum(b)
/* Weights matrix is now */
weights = b/c
st_matrix("weights2", weights)
end

* Most weight is given to the diarrhoea variable, followed by vomiting and then chills
foreach y of numlist 1/8 {
local we`y' = weights2[`y',1]
}

restore

replace morb_index = `we1'*std_sadl14_a + `we2'*std_carry_easy + `we3'*std_walk_easy + `we4'*std_sdiarrhoea14 + `we5'*std_sfever14 +`we6'*std_scough14 +`we7'*std_schills14 + `we8'*std_svomiting14 if male==0


global regressors2 "wg_o age edu_2 edu_3 edu_4 age2 last"

*-----------------------------------
*LINEAR PROBABILITY MODELS
*-------------------------------------
*Specify the smallest education level as the base level
fvset base 0 edu

matrix coeffm = J(1,9,.)
matrix colnames coeffm = "Health Index" "No diarrhoea" "No fever" "No cough" "No chills" "No vomiting" "Carry out daily activities" "Easily carry 10kg load for 20m" "Easily walk 5km"
matrix CI1 = J(2,9,.)
matrix colnames CI1 = "Health Index" "No diarrhoea" "No fever" "No cough" "No chills" "No vomiting" "Carry out daily activities" "Easily carry 10kg load for 20m" "Easily walk 5km"
matrix rownames CI1 = ll95 ul95 

matrix coefff = J(1,9,.)
matrix colnames coefff = "Health Index" "No diarrhoea" "No fever" "No cough" "No chills" "No vomiting" "Carry out daily activities" "Easily carry 10kg load for 20m" "Easily walk 5km"
matrix CIf = J(2,9,.)
matrix colnames CIf = "Health Index" "No diarrhoea" "No fever" "No cough" "No chills" "No vomiting" "Carry out daily activities" "Easily carry 10kg load for 20m" "Easily walk 5km"
matrix rownames CIf = ll95 ul95 

preserve
keep if male==1
qui regress $health
keep if e(sample)
local i=0
foreach y of varlist morb_index sdiarrhoea14 sfever14 scough14 schills14 svomiting14 sadl14_a carry_easy walk_easy {
	local ++ i
display "-------- Males ---------"
reg `y' $regressors2, cluster(zone_o)
matrix coeffm[1,`i']=e(b)[1,1]
boottest wg_o, boot(wild) seed(10101) bootcl(zone_o)
matrix CI1[1,`i']=r(CI)[1,1]
matrix CI1[2,`i']=r(CI)[1,2]


}
restore

preserve
keep if male==0
qui regress $health
keep if e(sample)
local i=0
foreach y of varlist morb_index sdiarrhoea14 sfever14 scough14 schills14 svomiting14 sadl14_a carry_easy walk_easy {
	local ++ i
reg `y' $regressors2, cluster(zone_o)
matrix coefff[1,`i']=e(b)[1,1]
boottest wg_o, boot(wild) seed(10101) bootcl(zone_o)
matrix CIf[1,`i']=r(CI)[1,1]
matrix CIf[2,`i']=r(CI)[1,2]
}

mat ren coeffm Males
mat ren coefff Females

coefplot (matrix(Males[1,]), ci(CI1)) || (matrix(Females[1,]), ci(CIf)),  xline(0)  ciopts(recast(rcap)) citop
graph save "$figures/adulthealth_treat_control", replace
graph export "$figures/adulthealth_treat_control.pdf", replace

log close
exit

